In [ ]:
%pylab notebook
import datetime as dt
import glob
import matplotlib.pyplot as plt
import matplotlib.dates as md
#from nose.tools import set_trace
from charistools.convertors import Dem2Hypsometry
from charistools.convertors import Modice2Hypsometry
from charistools.hypsometry import Hypsometry
from charistools.meltModels import TriSurfTempIndexMelt
from charistools.meltModels import ImshowTriSurfMelt
from charistools.meltModels import PlotTriSurfInput
from charistools.meltModels import PlotTriSurfMelt
from charistools.modelEnv import ModelEnv
from charistools.timeSeries import TimeSeries
import numpy as np
from netCDF4 import Dataset
import numpy as np
import pandas as pd
import rasterio
import re
import seaborn as sns
import os
sns.set()
sns.axes_style("darkgrid")
In [ ]:
configFile = '/Users/brodzik/ipython_notebooks/charis/calibration_modelEnv_config.ini'
myEnv = ModelEnv(tileConfigFile=configFile)
# reset basin mask directory and filename patterns slightly
print("Before:")
print(myEnv.tileConfig['input']['fixed']['basin_mask']['dir'])
print(myEnv.tileConfig['input']['fixed']['basin_mask']['pattern'])
print(myEnv.tileConfig['hypsometry']['area_by_elevation']['dir'])
print(myEnv.tileConfig['hypsometry']['modice_min05yr_by_elevation']['dir'])
myEnv.tileConfig['input']['fixed']['basin_mask']['dir'] = "%MODEL_TOP_DIR%/basins/major_basins_from_GRDC/MODIStiles"
myEnv.tileConfig['input']['fixed']['basin_mask']['pattern'] = "%DRAINAGEID%_%TILEID%.tif"
myEnv.tileConfig['hypsometry']['area_by_elevation']['dir'] = "%MODEL_TOP_DIR%/derived_hypsometries/REECv0/basin_areas"
myEnv.tileConfig['hypsometry']['modice_min05yr_by_elevation']['dir'] = "%MODEL_TOP_DIR%/derived_hypsometries/REECv0/basin_areas"
print("After:")
print(myEnv.tileConfig['input']['fixed']['basin_mask']['dir'])
print(myEnv.tileConfig['input']['fixed']['basin_mask']['pattern'])
print(myEnv.tileConfig['hypsometry']['area_by_elevation']['dir'])
print(myEnv.tileConfig['hypsometry']['modice_min05yr_by_elevation']['dir'])
In [ ]:
# Basins with endorheic holes filled in
#majorBasinIDs = ["SY_SyrDarya_at_TyumenAryk",
# "AM_AmuDarya_at_Chatly_noholes",
# "IN_Indus_at_Kotri_nolobe_gcs_noholes",
# "GA_Ganges_at_Paksey",
# "BR_Bramaputra_at_Bahadurabad_noholes"]
# GRDC Basins, Indus corrected for lobe
#majorBasinIDs = ["SY_SyrDarya_at_TyumenAryk",
# "AM_AmuDarya_at_Chatly",
# "IN_Indus_at_Kotri",
# "GA_Ganges_at_Paksey",
# "BR_Bramaputra_at_Bahadurabad"]
majorBasinIDs = [#"SY_SyrDarya_at_TyumenAryk",
#"AM_AmuDarya_at_Chatly",
#"AM_AmuDarya_at_Chatly_noholes",
#"IN_Indus_at_Kotri",
#"IN_Indus_at_Kotri_nolobe",
"IN_Indus_at_Kotri_nolobe_gcs_noholes"]
#"GA_Ganges_at_Paksey",
#"BR_Bramaputra_at_Bahadurabad",
#"BR_Bramaputra_at_Bahadurabad_noholes"]
majorBasinNames = [#"SY (GRDC)",
#"AM (GRDC)",
#"AM (GRDC holes filled)",
#"IN (GRDC)",
#"IN (GRDC no lobe)",
"IN (GRDC no lobe, holes filled)"]
#"GA (GRDC)",
#"BR (GRDC)",
#"BR (GRDC holes filled)"]
majorBasinIDs
In [ ]:
#print(ModelEnv.tileIDs_for_drainage(myEnv, drainageID=majorBasinIDs[4]))
#print(myEnv.hypsometry_filename(type='area_by_elevation',
# contour_m=50,
# drainageID=majorBasinIDs[0]))
#print(myEnv.hypsometry_filename(type='modice_min05yr_by_elevation',
# drainageID=majorBasinIDs[0], modice_nstrikes=3))
In [ ]:
verbose = False
contour_m = 100
In [ ]:
for id in majorBasinIDs:
outFile = myEnv.hypsometry_filename(type='area_by_elevation',
contour_m=contour_m,
drainageID=id)
print("DEM : %s" % outFile)
hyps = Dem2Hypsometry(id, myEnv, contour_m=contour_m,
outfile=outFile, decimal_places=4,
verbose=verbose)
print("Total masked area = %f" % hyps.data.sum().sum())
for n in np.arange(3):
outFile = myEnv.hypsometry_filename(type='modice_min05yr_by_elevation',
drainageID=id,
contour_m=contour_m,
modice_nstrikes=n+1)
print("MODICE: %d : %s" % (n, outFile))
hyps = Modice2Hypsometry(id, myEnv, contour_m=contour_m,
modice_nstrikes=n+1,
outfile=outFile,
decimal_places=4,
verbose=verbose)
In [ ]:
my_series = pd.Series(majorBasinIDs)
df = pd.DataFrame(my_series)
df.columns = ['Basin']
df.set_index('Basin', inplace=True)
df["Name"] = majorBasinNames
df
In [ ]:
outdir = "/Users/brodzik/projects/CHARIS/derived_hypsometries/REECv0/basin_areas"
In [ ]:
fig, axes = plt.subplots(1)
for id in majorBasinIDs:
outFile = myEnv.hypsometry_filename(type='area_by_elevation',
contour_m=contour_m,
drainageID=id)
print(outFile)
hyps = Hypsometry(outFile)
label = "%s (%d $km^2$)" % (id, hyps.data.sum().sum())
df.ix[id,'Basin Area (km2)'] = hyps.data.sum().sum()
axes.plot(hyps.data.loc['NoDate'].values,
np.array(hyps.data.loc['NoDate'].index),
label=label)
axes.legend()
fig.suptitle('CHARIS Major basins')
#fig.savefig("%s/%s" % (outdir, "CHARIS_basins_area_by_elevation.png"))
In [ ]:
df
In [ ]:
limit_m = 2000
fig, axes = plt.subplots()
for id in majorBasinIDs:
outFile = myEnv.hypsometry_filename(type='area_by_elevation',
contour_m = contour_m,
drainageID=id)
print(outFile)
hyps = Hypsometry(outFile)
hyps.data.drop(hyps.data.columns[hyps.data.columns < limit_m], axis=1, inplace=True)
label = "%s (%d $km^2$)" % (df.Name.loc[id], hyps.data.sum().sum())
df.ix[id,'Basin Area > %dm (km2)' % limit_m] = hyps.data.sum().sum()
axes.plot(hyps.data.loc['NoDate'].values,
np.array(hyps.data.loc['NoDate'].index),
label=label)
axes.set_xlabel("Area per %d m band ($km^2$)" % contour_m)
axes.set_ylabel("Elevation ($m$)")
axes.set_xlim([0., 25000.])
axes.set_ylim([2000., 7000.])
axes.legend()
fig.suptitle('Elevation distribution of CHARIS major basins')
#fig.savefig("%s/CHARIS_basin_area_by_elevation.above%dm.v1.png" % (outdir, limit_m))
In [ ]:
contour_m = 100
limit_m = 2000
for id in majorBasinIDs:
for n in np.arange(3) + 1:
outFile = myEnv.hypsometry_filename(type='modice_min05yr_by_elevation',
drainageID=id,
contour_m=contour_m,
modice_nstrikes=n)
print("MODICE: %d : %s" % (n, outFile))
hyps = Hypsometry(outFile)
print(hyps.data.transpose())
print(" Total MODICE area: %f" % hyps.data.sum().sum())
hyps.data.drop(hyps.data.columns[hyps.data.columns < limit_m],
axis=1, inplace=True)
print(" Total MODICE area (>2000): %f" % hyps.data.sum().sum())
df.ix[id,'MODICE (%dstr) Area > %dm (km2)' % (n, limit_m)] = hyps.data.sum().sum()
In [ ]:
df
In [ ]:
csvfile = "%s/CHARIS_basin_area_summary.csv" % outdir
df.to_csv(csvfile)
In [ ]:
outdir
In [ ]:
basin_dir = "/Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles"
modice_dir = "/Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc"
#drainageID = "SY_SyrDarya_at_TyumenAryk"
#tiles = ['h22v04', 'h23v04', 'h23v05']
#drainageID = "AM_AmuDarya_at_Chatly_noholes"
#tiles = ['h22v04', 'h22v05', 'h23v04', 'h23v05']
#drainageID = "IN_Indus_at_Kotri_nolobe_gcs_noholes"
#tiles = ['h23v05', 'h23v06', 'h24v05', 'h24v06', 'h25v05']
#drainageID = "GA_Ganges_at_Paksey"
#tiles = ['h24v05', 'h24v06', 'h25v05', 'h25v06', 'h26v06']
drainageID = "BR_Bramaputra_at_Bahadurabad_noholes"
tiles = ['h25v05', 'h25v06', 'h26v05', 'h26v06']
nstrikes = 2
In [ ]:
total_area = 0.
for tile in tiles:
basinFile = "%s/%s_%s.tif" % (basin_dir, drainageID, tile)
modiceFile = "%s/MODICE.v0.4.%s.%dstrike.min05yr.mask.nc" % (modice_dir, tile, nstrikes)
# Read mask file
with rasterio.open(basinFile) as src:
mask = np.squeeze(src.read())
print(basinFile)
f = Dataset(modiceFile, "r", 'NETCDF4')
modice = f.variables['modice_min_year_mask'][:]
modice[modice == 1] = 0
print(modiceFile)
modice_count = modice[modice == 2].shape[0]
unmasked_area = modice_count * 463.312717 * 463.312717 / ( 1000. * 1000.)
fig, axes = plt.subplots(2, 2, figsize=(12,8))
axes[0,0].imshow(mask, cmap='Greys', vmin=np.amin(mask), vmax=np.amax(mask), interpolation='None')
axes[0,0].set_title("%s: mask" % drainageID)
axes[0,0].axis('off')
axes[0,1].imshow(modice, cmap='Greys', vmin=np.amin(modice), vmax=np.amax(modice), interpolation='None')
axes[0,1].set_title("%s: modice" % tile )
axes[0,1].axis('off')
masked_modice = modice.copy()
masked_modice[mask == 0] = 0
masked_count = masked_modice[masked_modice == 2].shape[0]
masked_area = masked_count * 463.312717 * 463.312717 / ( 1000. * 1000.)
print("in %s, %s, ice pixels %d, area = %f sq km" % (drainageID, tile, masked_count, masked_area ))
axes[1,1].imshow(masked_modice, cmap='Greys', vmin=np.amin(modice), vmax=np.amax(modice), interpolation='None')
axes[1,1].set_title("%s: modice masked for %s" % (tile, drainageID) )
axes[1,1].axis('off')
total_area = total_area + masked_area
fig.tight_layout()
print("in %s, MODICE area = %f" % (drainageID, total_area))
In [ ]: